<!DOCTYPE html>
<HTML lang = "en">
<HEAD>
  <meta charset="UTF-8"/>
  <meta name="viewport" content="width=device-width, initial-scale=1.0, user-scalable=yes">
  <title>Global Sensitivity Analysis</title>
  

  <script type="text/x-mathjax-config">
    MathJax.Hub.Config({
      tex2jax: {inlineMath: [['$','$'], ['\\(','\\)']]},
      TeX: { equationNumbers: { autoNumber: "AMS" } }
    });
  </script>

  <script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
  </script>

  
<style>
pre.hljl {
    border: 1px solid #ccc;
    margin: 5px;
    padding: 5px;
    overflow-x: auto;
    color: rgb(68,68,68); background-color: rgb(251,251,251); }
pre.hljl > span.hljl-t { }
pre.hljl > span.hljl-w { }
pre.hljl > span.hljl-e { }
pre.hljl > span.hljl-eB { }
pre.hljl > span.hljl-o { }
pre.hljl > span.hljl-k { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kc { color: rgb(59,151,46); font-style: italic; }
pre.hljl > span.hljl-kd { color: rgb(214,102,97); font-style: italic; }
pre.hljl > span.hljl-kn { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kp { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kr { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-kt { color: rgb(148,91,176); font-weight: bold; }
pre.hljl > span.hljl-n { }
pre.hljl > span.hljl-na { }
pre.hljl > span.hljl-nb { }
pre.hljl > span.hljl-nbp { }
pre.hljl > span.hljl-nc { }
pre.hljl > span.hljl-ncB { }
pre.hljl > span.hljl-nd { color: rgb(214,102,97); }
pre.hljl > span.hljl-ne { }
pre.hljl > span.hljl-neB { }
pre.hljl > span.hljl-nf { color: rgb(66,102,213); }
pre.hljl > span.hljl-nfm { color: rgb(66,102,213); }
pre.hljl > span.hljl-np { }
pre.hljl > span.hljl-nl { }
pre.hljl > span.hljl-nn { }
pre.hljl > span.hljl-no { }
pre.hljl > span.hljl-nt { }
pre.hljl > span.hljl-nv { }
pre.hljl > span.hljl-nvc { }
pre.hljl > span.hljl-nvg { }
pre.hljl > span.hljl-nvi { }
pre.hljl > span.hljl-nvm { }
pre.hljl > span.hljl-l { }
pre.hljl > span.hljl-ld { color: rgb(148,91,176); font-style: italic; }
pre.hljl > span.hljl-s { color: rgb(201,61,57); }
pre.hljl > span.hljl-sa { color: rgb(201,61,57); }
pre.hljl > span.hljl-sb { color: rgb(201,61,57); }
pre.hljl > span.hljl-sc { color: rgb(201,61,57); }
pre.hljl > span.hljl-sd { color: rgb(201,61,57); }
pre.hljl > span.hljl-sdB { color: rgb(201,61,57); }
pre.hljl > span.hljl-sdC { color: rgb(201,61,57); }
pre.hljl > span.hljl-se { color: rgb(59,151,46); }
pre.hljl > span.hljl-sh { color: rgb(201,61,57); }
pre.hljl > span.hljl-si { }
pre.hljl > span.hljl-so { color: rgb(201,61,57); }
pre.hljl > span.hljl-sr { color: rgb(201,61,57); }
pre.hljl > span.hljl-ss { color: rgb(201,61,57); }
pre.hljl > span.hljl-ssB { color: rgb(201,61,57); }
pre.hljl > span.hljl-nB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nbB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nfB { color: rgb(59,151,46); }
pre.hljl > span.hljl-nh { color: rgb(59,151,46); }
pre.hljl > span.hljl-ni { color: rgb(59,151,46); }
pre.hljl > span.hljl-nil { color: rgb(59,151,46); }
pre.hljl > span.hljl-noB { color: rgb(59,151,46); }
pre.hljl > span.hljl-oB { color: rgb(102,102,102); font-weight: bold; }
pre.hljl > span.hljl-ow { color: rgb(102,102,102); font-weight: bold; }
pre.hljl > span.hljl-p { }
pre.hljl > span.hljl-c { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-ch { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cm { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cp { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cpB { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-cs { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-csB { color: rgb(153,153,119); font-style: italic; }
pre.hljl > span.hljl-g { }
pre.hljl > span.hljl-gd { }
pre.hljl > span.hljl-ge { }
pre.hljl > span.hljl-geB { }
pre.hljl > span.hljl-gh { }
pre.hljl > span.hljl-gi { }
pre.hljl > span.hljl-go { }
pre.hljl > span.hljl-gp { }
pre.hljl > span.hljl-gs { }
pre.hljl > span.hljl-gsB { }
pre.hljl > span.hljl-gt { }
</style>



  <style type="text/css">
  @font-face {
  font-style: normal;
  font-weight: 300;
}
@font-face {
  font-style: normal;
  font-weight: 400;
}
@font-face {
  font-style: normal;
  font-weight: 600;
}
html {
  font-family: sans-serif; /* 1 */
  -ms-text-size-adjust: 100%; /* 2 */
  -webkit-text-size-adjust: 100%; /* 2 */
}
body {
  margin: 0;
}
article,
aside,
details,
figcaption,
figure,
footer,
header,
hgroup,
main,
menu,
nav,
section,
summary {
  display: block;
}
audio,
canvas,
progress,
video {
  display: inline-block; /* 1 */
  vertical-align: baseline; /* 2 */
}
audio:not([controls]) {
  display: none;
  height: 0;
}
[hidden],
template {
  display: none;
}
a:active,
a:hover {
  outline: 0;
}
abbr[title] {
  border-bottom: 1px dotted;
}
b,
strong {
  font-weight: bold;
}
dfn {
  font-style: italic;
}
h1 {
  font-size: 2em;
  margin: 0.67em 0;
}
mark {
  background: #ff0;
  color: #000;
}
small {
  font-size: 80%;
}
sub,
sup {
  font-size: 75%;
  line-height: 0;
  position: relative;
  vertical-align: baseline;
}
sup {
  top: -0.5em;
}
sub {
  bottom: -0.25em;
}
img {
  border: 0;
}
svg:not(:root) {
  overflow: hidden;
}
figure {
  margin: 1em 40px;
}
hr {
  -moz-box-sizing: content-box;
  box-sizing: content-box;
  height: 0;
}
pre {
  overflow: auto;
}
code,
kbd,
pre,
samp {
  font-family: monospace, monospace;
  font-size: 1em;
}
button,
input,
optgroup,
select,
textarea {
  color: inherit; /* 1 */
  font: inherit; /* 2 */
  margin: 0; /* 3 */
}
button {
  overflow: visible;
}
button,
select {
  text-transform: none;
}
button,
html input[type="button"], /* 1 */
input[type="reset"],
input[type="submit"] {
  -webkit-appearance: button; /* 2 */
  cursor: pointer; /* 3 */
}
button[disabled],
html input[disabled] {
  cursor: default;
}
button::-moz-focus-inner,
input::-moz-focus-inner {
  border: 0;
  padding: 0;
}
input {
  line-height: normal;
}
input[type="checkbox"],
input[type="radio"] {
  box-sizing: border-box; /* 1 */
  padding: 0; /* 2 */
}
input[type="number"]::-webkit-inner-spin-button,
input[type="number"]::-webkit-outer-spin-button {
  height: auto;
}
input[type="search"] {
  -webkit-appearance: textfield; /* 1 */
  -moz-box-sizing: content-box;
  -webkit-box-sizing: content-box; /* 2 */
  box-sizing: content-box;
}
input[type="search"]::-webkit-search-cancel-button,
input[type="search"]::-webkit-search-decoration {
  -webkit-appearance: none;
}
fieldset {
  border: 1px solid #c0c0c0;
  margin: 0 2px;
  padding: 0.35em 0.625em 0.75em;
}
legend {
  border: 0; /* 1 */
  padding: 0; /* 2 */
}
textarea {
  overflow: auto;
}
optgroup {
  font-weight: bold;
}
table {
  font-family: monospace, monospace;
  font-size : 0.8em;
  border-collapse: collapse;
  border-spacing: 0;
}
td,
th {
  padding: 0;
}
thead th {
    border-bottom: 1px solid black;
    background-color: white;
}
tr:nth-child(odd){
  background-color: rgb(248,248,248);
}


/*
* Skeleton V2.0.4
* Copyright 2014, Dave Gamache
* www.getskeleton.com
* Free to use under the MIT license.
* http://www.opensource.org/licenses/mit-license.php
* 12/29/2014
*/
.container {
  position: relative;
  width: 100%;
  max-width: 960px;
  margin: 0 auto;
  padding: 0 20px;
  box-sizing: border-box; }
.column,
.columns {
  width: 100%;
  float: left;
  box-sizing: border-box; }
@media (min-width: 400px) {
  .container {
    width: 85%;
    padding: 0; }
}
@media (min-width: 550px) {
  .container {
    width: 80%; }
  .column,
  .columns {
    margin-left: 4%; }
  .column:first-child,
  .columns:first-child {
    margin-left: 0; }

  .one.column,
  .one.columns                    { width: 4.66666666667%; }
  .two.columns                    { width: 13.3333333333%; }
  .three.columns                  { width: 22%;            }
  .four.columns                   { width: 30.6666666667%; }
  .five.columns                   { width: 39.3333333333%; }
  .six.columns                    { width: 48%;            }
  .seven.columns                  { width: 56.6666666667%; }
  .eight.columns                  { width: 65.3333333333%; }
  .nine.columns                   { width: 74.0%;          }
  .ten.columns                    { width: 82.6666666667%; }
  .eleven.columns                 { width: 91.3333333333%; }
  .twelve.columns                 { width: 100%; margin-left: 0; }

  .one-third.column               { width: 30.6666666667%; }
  .two-thirds.column              { width: 65.3333333333%; }

  .one-half.column                { width: 48%; }

  /* Offsets */
  .offset-by-one.column,
  .offset-by-one.columns          { margin-left: 8.66666666667%; }
  .offset-by-two.column,
  .offset-by-two.columns          { margin-left: 17.3333333333%; }
  .offset-by-three.column,
  .offset-by-three.columns        { margin-left: 26%;            }
  .offset-by-four.column,
  .offset-by-four.columns         { margin-left: 34.6666666667%; }
  .offset-by-five.column,
  .offset-by-five.columns         { margin-left: 43.3333333333%; }
  .offset-by-six.column,
  .offset-by-six.columns          { margin-left: 52%;            }
  .offset-by-seven.column,
  .offset-by-seven.columns        { margin-left: 60.6666666667%; }
  .offset-by-eight.column,
  .offset-by-eight.columns        { margin-left: 69.3333333333%; }
  .offset-by-nine.column,
  .offset-by-nine.columns         { margin-left: 78.0%;          }
  .offset-by-ten.column,
  .offset-by-ten.columns          { margin-left: 86.6666666667%; }
  .offset-by-eleven.column,
  .offset-by-eleven.columns       { margin-left: 95.3333333333%; }

  .offset-by-one-third.column,
  .offset-by-one-third.columns    { margin-left: 34.6666666667%; }
  .offset-by-two-thirds.column,
  .offset-by-two-thirds.columns   { margin-left: 69.3333333333%; }

  .offset-by-one-half.column,
  .offset-by-one-half.columns     { margin-left: 52%; }

}
html {
  font-size: 62.5%; }
body {
  font-size: 1.5em; /* currently ems cause chrome bug misinterpreting rems on body element */
  line-height: 1.6;
  font-weight: 400;
  font-family: "Raleway", "HelveticaNeue", "Helvetica Neue", Helvetica, Arial, sans-serif;
  color: #222; }
h1, h2, h3, h4, h5, h6 {
  margin-top: 0;
  margin-bottom: 2rem;
  font-weight: 300; }
h1 { font-size: 3.6rem; line-height: 1.2;  letter-spacing: -.1rem;}
h2 { font-size: 3.4rem; line-height: 1.25; letter-spacing: -.1rem; }
h3 { font-size: 3.2rem; line-height: 1.3;  letter-spacing: -.1rem; }
h4 { font-size: 2.8rem; line-height: 1.35; letter-spacing: -.08rem; }
h5 { font-size: 2.4rem; line-height: 1.5;  letter-spacing: -.05rem; }
h6 { font-size: 1.5rem; line-height: 1.6;  letter-spacing: 0; }

p {
  margin-top: 0; }
a {
  color: #1EAEDB; }
a:hover {
  color: #0FA0CE; }
.button,
button,
input[type="submit"],
input[type="reset"],
input[type="button"] {
  display: inline-block;
  height: 38px;
  padding: 0 30px;
  color: #555;
  text-align: center;
  font-size: 11px;
  font-weight: 600;
  line-height: 38px;
  letter-spacing: .1rem;
  text-transform: uppercase;
  text-decoration: none;
  white-space: nowrap;
  background-color: transparent;
  border-radius: 4px;
  border: 1px solid #bbb;
  cursor: pointer;
  box-sizing: border-box; }
.button:hover,
button:hover,
input[type="submit"]:hover,
input[type="reset"]:hover,
input[type="button"]:hover,
.button:focus,
button:focus,
input[type="submit"]:focus,
input[type="reset"]:focus,
input[type="button"]:focus {
  color: #333;
  border-color: #888;
  outline: 0; }
.button.button-primary,
button.button-primary,
input[type="submit"].button-primary,
input[type="reset"].button-primary,
input[type="button"].button-primary {
  color: #FFF;
  background-color: #33C3F0;
  border-color: #33C3F0; }
.button.button-primary:hover,
button.button-primary:hover,
input[type="submit"].button-primary:hover,
input[type="reset"].button-primary:hover,
input[type="button"].button-primary:hover,
.button.button-primary:focus,
button.button-primary:focus,
input[type="submit"].button-primary:focus,
input[type="reset"].button-primary:focus,
input[type="button"].button-primary:focus {
  color: #FFF;
  background-color: #1EAEDB;
  border-color: #1EAEDB; }
input[type="email"],
input[type="number"],
input[type="search"],
input[type="text"],
input[type="tel"],
input[type="url"],
input[type="password"],
textarea,
select {
  height: 38px;
  padding: 6px 10px; /* The 6px vertically centers text on FF, ignored by Webkit */
  background-color: #fff;
  border: 1px solid #D1D1D1;
  border-radius: 4px;
  box-shadow: none;
  box-sizing: border-box; }
/* Removes awkward default styles on some inputs for iOS */
input[type="email"],
input[type="number"],
input[type="search"],
input[type="text"],
input[type="tel"],
input[type="url"],
input[type="password"],
textarea {
  -webkit-appearance: none;
     -moz-appearance: none;
          appearance: none; }
textarea {
  min-height: 65px;
  padding-top: 6px;
  padding-bottom: 6px; }
input[type="email"]:focus,
input[type="number"]:focus,
input[type="search"]:focus,
input[type="text"]:focus,
input[type="tel"]:focus,
input[type="url"]:focus,
input[type="password"]:focus,
textarea:focus,
select:focus {
  border: 1px solid #33C3F0;
  outline: 0; }
label,
legend {
  display: block;
  margin-bottom: .5rem;
  font-weight: 600; }
fieldset {
  padding: 0;
  border-width: 0; }
input[type="checkbox"],
input[type="radio"] {
  display: inline; }
label > .label-body {
  display: inline-block;
  margin-left: .5rem;
  font-weight: normal; }
ul {
  list-style: circle; }
ol {
  list-style: decimal; }
ul ul,
ul ol,
ol ol,
ol ul {
  margin: 1.5rem 0 1.5rem 3rem;
  font-size: 90%; }
li > p {margin : 0;}
th,
td {
  padding: 12px 15px;
  text-align: left;
  border-bottom: 1px solid #E1E1E1; }
th:first-child,
td:first-child {
  padding-left: 0; }
th:last-child,
td:last-child {
  padding-right: 0; }
button,
.button {
  margin-bottom: 1rem; }
input,
textarea,
select,
fieldset {
  margin-bottom: 1.5rem; }
pre,
blockquote,
dl,
figure,
table,
p,
ul,
ol,
form {
  margin-bottom: 1.0rem; }
.u-full-width {
  width: 100%;
  box-sizing: border-box; }
.u-max-full-width {
  max-width: 100%;
  box-sizing: border-box; }
.u-pull-right {
  float: right; }
.u-pull-left {
  float: left; }
hr {
  margin-top: 3rem;
  margin-bottom: 3.5rem;
  border-width: 0;
  border-top: 1px solid #E1E1E1; }
.container:after,
.row:after,
.u-cf {
  content: "";
  display: table;
  clear: both; }

pre {
  display: block;
  padding: 9.5px;
  margin: 0 0 10px;
  font-size: 13px;
  line-height: 1.42857143;
  word-break: break-all;
  word-wrap: break-word;
  border: 1px solid #ccc;
  border-radius: 4px;
}

pre.hljl {
  margin: 0 0 10px;
  display: block;
  background: #f5f5f5;
  border-radius: 4px;
  padding : 5px;
}

pre.output {
  background: #ffffff;
}

pre.code {
  background: #ffffff;
}

pre.julia-error {
  color : red
}

code,
kbd,
pre,
samp {
  font-family: Menlo, Monaco, Consolas, "Courier New", monospace;
  font-size: 0.9em;
}


@media (min-width: 400px) {}
@media (min-width: 550px) {}
@media (min-width: 750px) {}
@media (min-width: 1000px) {}
@media (min-width: 1200px) {}

h1.title {margin-top : 20px}
img {max-width : 100%}
div.title {text-align: center;}

  </style>
</HEAD>

<BODY>
  <div class ="container">
    <div class = "row">
      <div class = "col-md-12 twelve columns">
        <div class="title">
          <h1 class="title">Global Sensitivity Analysis</h1>
          <h5>Chris Rackauckas</h5>
          <h5>December 12st, 2020</h5>
        </div>

        <p>Sensitivity analysis is the measure of how sensitive a model is to changes in parameters, i.e. how much the output changes given a change in the input. Clearly, derivatives are a measure of sensitivity, but derivative are <em>local sensitivity</em> measures because they are only the derivative at a single point. However, the idea of probabilistic programming starts to bring up an alternative question: how does the output of a model generally change with a change in the input? This kind of question requires an understanding of <em>global sensitivty</em> of a model. While there isn&#39;t a single definition of the concept, there are a few methods that individuals have employed to estimate the global sensitivity.</p>
<p>Reference implementations of these methods can be found in <a href="https://github.com/SciML/GlobalSensitivity.jl">GlobalSensitivity.jl</a></p>
<h2>Setup for Global Sensitivity</h2>
<p>In our global sensitivity analysis, we have a model <span class="math">$f$</span> and want to understand the relationship</p>
<p class="math">\[
y = f(x_i)
\]</p>
<p>Recall <span class="math">$f$</span> can be a neural network, an ODE solve, etc. where the <span class="math">$X_i$</span> are items like initial conditions and parameters. What we want to do is understand how much the total changes in <span class="math">$y$</span> can be attributed to changes in specific <span class="math">$x_i$</span>.</p>
<p>However, this is not an actionable form since we don&#39;t know what valid inputs into <span class="math">$f$</span> look like. Thus any global sensitivity study at least needs a domain for the <span class="math">$x_i$</span>, at least in terms of bounds. This is still underdefined because what makes one thing that it&#39;s not more likely for <span class="math">$x_i$</span> to be near the lower part of the bound instead of the upper part? Thus, for global sensitivity analysis to be well-defined, <span class="math">$x_i$</span> must take a distributional form, i.e. be random variables. Thus <span class="math">$f$</span> is a determinstic program with probabilistic inputs, and we want to determine the effects of the distributional inputs on the distribution of the output.</p>
<h3>Reasons for Global Sensitivity Analysis</h3>
<p>What are the things we can learn from doing such a global sensitivity analysis?</p>
<ol>
<li><p>You can learn what variables would need to be changed to drive the solution in a given direction or control the system. If your model is exact and the parameters are known, the &quot;standard&quot; methods apply, but if your model is only approximate, a global sensitivity metric may be a better prediction as to how variables cause changes.</p>
</li>
<li><p>You can learn if there are any variables which do not have a true effect on the output. These variables would be practically unidentifiable from data and models can be reduced by removing the terms. It also is predictive as to robustness properties.</p>
</li>
<li><p>You can find ways to automatically sparsify a model by dropping off the components which contribute the least. This matters in automatically generated or automatically detected models, where many pieces may be spurious and global sensitivities would be a method to detect that in a manner that is not sensitive to the chosen parameters.</p>
</li>
</ol>
<h2>Global Sensitivity Analysis Measures</h2>
<h3>Linear Global Sensitivity Metrics: Correlations and Regressions</h3>
<p>The first thing that you can do is approximate the full model with a linear surrogate, i.e.</p>
<p class="math">\[
y = AX
\]</p>
<p>for some linear model. A regression can be done on the outputs of the model in order to find the linear approximation. The best fitting global linear model then gives coefficients for the global sensitivities via the individual effects, i.e. for</p>
<p class="math">\[
y = \sum_i \beta_i x_i
\]</p>
<p>,</p>
<p>the <span class="math">$\beta_i$</span> are the global effect. Just as with any use of a linear model, the same ideas apply. The coefficient of determination &#40;<span class="math">$R^2$</span>&#41; is a measure of how well the model fits. However, one major change needs to be done in order to ensure that the solutions are comparable between different models. The dependence of the solution on the units can cause the coefficients to be large/small. Thus we need to normalize the data, i.e. use the transformation</p>
<p class="math">\[
\tilde{x_i} = \frac{x_i-E[x_i]}{V[x_i]}
\]</p>
<p class="math">\[
\tilde{y_i} = \frac{y_i-E[y_i]}{V[y_i]}
\]</p>
<p>The normalized coefficients are known as the <em>Standardized Regression Coefficients</em> &#40;SRC&#41; and are a measure of the global effects.</p>
<p>Notice that while the <span class="math">$\beta_i$</span> capture the mean effects, it holds that</p>
<p class="math">\[
V(y) = \sum_i \beta^2_i x_i
\]</p>
<p>and thus the variance due to <span class="math">$x_i$</span> can be measured as:</p>
<p class="math">\[
SRC_i = \beta_i \sqrt{\frac{V[x_i]}{V[y]}}
\]</p>
<p>This interpretation is the same as the solution from the normalized variables.</p>
<p>From the same linear model, two other global sensitivity metrics are defined. The <em>Correlation Coefficients</em> &#40;CC&#41; are simply the correlations:</p>
<p class="math">\[
CC_i = \frac{\text{cov}(x_i,y)}{\sqrt{V[x_i]V[y]}}
\]</p>
<p>Similarly, the <em>Partial Correlation Coefficient</em> is the correlation coefficient where the linear effect of the other terms are removed, i.e. for <span class="math">$S_i = {x_1,x_2,\ldots,x_{j-1},x_{j+1},\ldots,x_n}$</span> we have</p>
<p class="math">\[
PCC_{i|S_i} = \frac{\text{cov}(x_i,y|S)j)}{\sqrt{V[x_i|S_i]V[y|S_i]}}
\]</p>
<h3>Derivative-based Global Sensitivity Measures &#40;DGSM&#41;</h3>
<p>To go beyond just a linear model, one might want to do successive linearization. Since derivatives are a form of linearization, then one may thing to average derivatives. This averaging of derivatives is the DGSM method. If the <span class="math">$x_i$</span> are random variables with joint CDF <span class="math">$F(x)$</span>, then it holds that:</p>
<p class="math">\[
v_i = \int_{R^d} \left(\frac{\partial f(x)}{\partial x_i}\right)^2 dF(x) = \mathbb{E}\left[\left(\frac{\partial f(x)}{\partial x_i}\right)^2\right],
\]</p>
<p>We can also define the mean measure, which is simply:</p>
<p class="math">\[
w_i = \int_{R^d} \frac{\partial f(x)}{\partial x_i} dF(x) = \mathbb{E}\left[\frac{\partial f(x)}{\partial x_i}\right].
\]</p>
<p>Thus a global variance estimate would be <span class="math">$v_i - w_i^2$</span>.</p>
<h3>ADVI for Global Sensitivity</h3>
<p>Note that the previously discussed method for probabilistic programming, ADVI, is a method for producing a Gaussian approximation for a probabilistic program. The resulting mean-field or full Gaussian approximations are variance index calculations&#33;</p>
<h3>The Morris One-At-A-Time &#40;OAT&#41; Method</h3>
<p>Instead of using derivatives, one can use finite difference approximations. Normally you want to use small <span class="math">$\Delta x$</span>, but if we are averaging derivatives over a large area, then in reality we don&#39;t really need a small <span class="math">$\Delta x$</span>&#33;</p>
<p>This is where the Morris method comes in. The basic idea is that moving in one direction at a time is a derivative estimate, and if we step large enough then the next derivative estimate may be sufficiently different enough to contribute well to the total approximation. Thus we do the following:</p>
<ol>
<li><p>Take a random starting point</p>
</li>
<li><p>Randomly choose a direction <span class="math">$i$</span> and make a change <span class="math">$\Delta x_i$</span> only in that direction.</p>
</li>
<li><p>Calculate the derivative approximation from that change. Repeat 2 and 3.</p>
</li>
</ol>
<p>Keep doing this for enough steps, and the average of your derivative approximations becomes a global index. Notice that this reuses every simulation as part of two separate estimates, making it much more computationally efficient than the other methods. However, it accounts for average changes and not necessarily measurements gives a value that&#39;s a decomposition of a total variance. But its computational cost makes it attractive for making quick estimates of the global sensitivities.</p>
<p>For practical usage, a few changes have to be done. First of all, notice that positive and negative change can cancel out. Thus if one wanta measure of associated variance, one should use absolute values or squared differences. Also, one needs to make sure that these trajectories get good coverage of the input space. Define the distance between two trajectories as the sum of the geometric distances between all pairs of points. Generate many more trajectories than necessary and choose the <span class="math">$r$</span> trajectories with the largest distance. If the model evaluations are expensive, this is significantly cheap enough in comparison that it&#39;s okay to do.</p>
<h3>Sobol&#39;s Method &#40;ANOVA&#41;</h3>
<p>Sobol&#39;s method is a true nonlinear decomposition of variance and it is thus considered one of the gold standards. For Sobol&#39;s method, we define the decomposition</p>
<p class="math">\[
f(x) = f_0 + \sum_i f_i(x_i) + \sum_{i,j} f_{ij}(x_i,x_j) + \ldots
\]</p>
<p>where</p>
<p class="math">\[
f_0 = \int_\Omega f(x) dx
\]</p>
<p>and orthogonality holds:</p>
<p class="math">\[
f_{i,j,\ldots}(x_i,x_j,\ldots)dx = 0
\]</p>
<p>by the definitions:</p>
<p class="math">\[
f_i(x_i) = E(y|x_i) - f_0
\]</p>
<p class="math">\[
f_{ij}(x_i,y_j) = E(y|x_i,x_j) - f_0 - f_i - f_j
\]</p>
<p>Assuming that <span class="math">$f(x)$</span> is L2, it holds that</p>
<p>\int<em>\Omega f^2&#40;x&#41;dx - f</em>0^2 &#61; \sum<em>s \sum</em>i \int f^2<em>&#123;i</em>1,i<em>2,\ldots,i</em>s&#125; dx&#36;</p>
<p>and thus</p>
<p class="math">\[
V[y] = \sum V_i + \sum V_{ij} + \ldots
\]</p>
<p>where</p>
<p class="math">\[
V_i = V[E_{x_{\sim i}}[y|x_i]]
\]</p>
<p class="math">\[
V_{ij} = V[E_{x_{\sim ij}}[y|x_i,x_j]]-V_i - V_j
\]</p>
<p>where <span class="math">$X_{\sim i}$</span> means all of the variables except <span class="math">$X_i$</span>. This means that the total variance can be decomposed into each of these variances.</p>
<p>From there, the fractional contribution to the total variance is thus the index:</p>
<p class="math">\[
S_i = \frac{V_i}{Var[y]}
\]</p>
<p>and similary for the second, third, etc. indices.</p>
<p>Additionally, if there are too many variables, once can compute the contribution of <span class="math">$x_i$</span> including all of its interactions as:</p>
<p class="math">\[
S_{T_i} = \frac{E_{X_{\sim i}}[Var[y|X_{\sim i}]]}{Var[y]} = 1 - \frac{Var_{X_{\sim i}}[E_{X_i}[y|x_{\sim i}]]}{Var[y]}
\]</p>
<h4>Computational Tractability and Quasi-Monte Carlo</h4>
<p>Notice that every single expectation has an integral in it, so the variance is defined as integrals of integrals, making this a very challenging calculation. Thus instead of directly calculating the integrals, in many cases Monte Carlo estimators are used. Instead of a pure Monte Carlo method, one generally uses a low-discrepency sequence &#40;a form of quasi-Monte Carlo&#41; to effectively sample the search space.</p>
<p>The following generates for example a <em>Sobol sequence</em>:</p>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>Sobol</span><span class='hljl-p'>,</span><span class='hljl-t'> </span><span class='hljl-n'>Plots</span><span class='hljl-t'>
</span><span class='hljl-n'>s</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>SobolSeq</span><span class='hljl-p'>(</span><span class='hljl-ni'>2</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-n'>p</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>hcat</span><span class='hljl-p'>([</span><span class='hljl-nf'>next!</span><span class='hljl-p'>(</span><span class='hljl-n'>s</span><span class='hljl-p'>)</span><span class='hljl-t'> </span><span class='hljl-k'>for</span><span class='hljl-t'> </span><span class='hljl-n'>i</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-ni'>1</span><span class='hljl-oB'>:</span><span class='hljl-ni'>1024</span><span class='hljl-p'>]</span><span class='hljl-oB'>...</span><span class='hljl-p'>)</span><span class='hljl-oB'>&#39;</span><span class='hljl-t'>
</span><span class='hljl-nf'>scatter</span><span class='hljl-p'>(</span><span class='hljl-n'>p</span><span class='hljl-p'>[</span><span class='hljl-oB'>:</span><span class='hljl-p'>,</span><span class='hljl-ni'>1</span><span class='hljl-p'>],</span><span class='hljl-t'> </span><span class='hljl-n'>p</span><span class='hljl-p'>[</span><span class='hljl-oB'>:</span><span class='hljl-p'>,</span><span class='hljl-ni'>2</span><span class='hljl-p'>])</span>
</pre>


<img src=""  />

<p>Another common quasi-Monte Carlo sequence is the <em>Latin Hypercube</em>, which is a generalization of the Latin Square where in every row, column, etc. only one point is given, allowing a linear spread over a high dimensional space.</p>


<pre class='hljl'>
<span class='hljl-k'>using</span><span class='hljl-t'> </span><span class='hljl-n'>LatinHypercubeSampling</span><span class='hljl-t'>
</span><span class='hljl-n'>p</span><span class='hljl-t'> </span><span class='hljl-oB'>=</span><span class='hljl-t'> </span><span class='hljl-nf'>LHCoptim</span><span class='hljl-p'>(</span><span class='hljl-ni'>120</span><span class='hljl-p'>,</span><span class='hljl-ni'>2</span><span class='hljl-p'>,</span><span class='hljl-ni'>1000</span><span class='hljl-p'>)</span><span class='hljl-t'>
</span><span class='hljl-nf'>scatter</span><span class='hljl-p'>(</span><span class='hljl-n'>p</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>][</span><span class='hljl-oB'>:</span><span class='hljl-p'>,</span><span class='hljl-ni'>1</span><span class='hljl-p'>],</span><span class='hljl-n'>p</span><span class='hljl-p'>[</span><span class='hljl-ni'>1</span><span class='hljl-p'>][</span><span class='hljl-oB'>:</span><span class='hljl-p'>,</span><span class='hljl-ni'>2</span><span class='hljl-p'>])</span>
</pre>


<img src=""  />

<p>For a reference library with many different quasi-Monte Carlo samplers, check out <a href="https://github.com/SciML/QuasiMonteCarlo.jl">QuasiMonteCarlo.jl</a>.</p>
<h2>Fourier Amplitude Sensitivity Sampling &#40;FAST&#41; and eFAST</h2>
<p>The FAST method is a change to the Sobol method to allow for faster convergence. First transform the variables <span class="math">$x_i$</span> onto the space <span class="math">$[0,1]$</span>. Then, instead of the linear decomposition, one decomposes into a Fourier basis:</p>
<p class="math">\[
f(x_i,x_2,\ldots,x_n) = \sum_{m_1 = -\infty}^{\infty} \ldots \sum_{m_n = -\infty}^{\infty} C_{m_1m_2\ldots m_n}\exp\left(2\pi i (m_1 x_1 + \ldots + m_n x_n)\right)
\]</p>
<p>where</p>
<p class="math">\[
C_{m_1m_2\ldots m_n} = \int_0^1 \ldots \int_0^1 f(x_i,x_2,\ldots,x_n) \exp\left(-2\pi i (m_1 x_1 + \ldots + m_n x_n)\right)
\]</p>
<p>The ANOVA like decomposition is thus</p>
<p class="math">\[
f_0 = C_{0\ldots 0}
\]</p>
<p class="math">\[
f_j = \sum_{m_j \neq 0} C_{0\ldots 0 m_j 0 \ldots 0} \exp (2\pi i m_j x_j)
\]</p>
<p class="math">\[
f_{jk} = \sum_{m_j \neq 0} \sum_{m_k \neq 0} C_{0\ldots 0 m_j 0 \ldots m_k 0 \ldots 0} \exp \left(2\pi i (m_j x_j + m_k x_k)\right)
\]</p>
<p>The first order conditional variance is thus:</p>
<p class="math">\[
V_j = \int_0^1 f_j^2 (x_j) dx_j = \sum_{m_j \neq 0} |C_{0\ldots 0 m_j 0 \ldots 0}|^2
\]</p>
<p>or</p>
<p class="math">\[
V_j = 2\sum_{m_j = 1}^\infty \left(A_{m_j}^2 + B_{m_j}^2 \right)
\]</p>
<p>where <span class="math">$C_{0\ldots 0 m_j 0 \ldots 0} = A_{m_j} + i B_{m_j}$</span>. By Fourier series we know this to be:</p>
<p class="math">\[
A_{m_j} = \int_0^1 \ldots \int_0^1 f(x)\cos(2\pi m_j x_j)dx
\]</p>
<p class="math">\[
B_{m_j} = \int_0^1 \ldots \int_0^1 f(x)\sin(2\pi m_j x_j)dx
\]</p>
<h4>Implementation via the Ergodic Theorem</h4>
<p>Define</p>
<p class="math">\[
X_j(s) = \frac{1}{2\pi} (\omega_j s \mod 2\pi)
\]</p>
<p>By the ergodic theorem, if <span class="math">$\omega_j$</span> are irrational numbers, then the dynamical system will never repeat values and thus it will create a solution that is dense in the plane &#40;Let&#39;s prove a bit later&#41;. As an animation:</p>
<p><img src="https://upload.wikimedia.org/wikipedia/commons/6/64/Search_curve_1.gif" alt="" /></p>
<p>&#40;here, <span class="math">$\omega_1 = \pi$</span> and <span class="math">$\omega_2 = 7$</span>&#41;</p>
<p>This means that:</p>
<p class="math">\[
A_{m_j} = \lim_{T\rightarrow \infty} \frac{1}{2T} \int_{-T}^T f(x)\cos(m_j \omega_j s)ds
\]</p>
<p class="math">\[
B_{m_j} = \lim_{T\rightarrow \infty} \frac{1}{2T} \int_{-T}^T f(x)\sin(m_j \omega_j s)ds
\]</p>
<p>i.e. the multidimensional integral can be approximated by the integral over a single line.</p>
<p>One can satisfy this approximately to get a simpler form for the integral. Using <span class="math">$\omega_i$</span> as integers, the integral is periodic and so only integrating over <span class="math">$2\pi$</span> is required. This would mean that:</p>
<p class="math">\[
A_{m_j} \approx \frac{1}{2\pi} \int_{-\pi}^\pi f(x)\cos(m_j \omega_j s)ds
\]</p>
<p class="math">\[
B_{m_j} \approx \frac{1}{2\pi} \int_{-\pi}^\pi f(x)\sin(m_j \omega_j s)ds
\]</p>
<p>It&#39;s only approximate since the sequence cannot be dense. For example, with <span class="math">$\omega_1 = 11$</span> and <span class="math">$\omega_2 = 7$</span>:</p>
<p><img src="https://upload.wikimedia.org/wikipedia/commons/2/29/Search_curve_3.gif" alt="" /></p>
<p>A higher period thus gives a better fill of the space and thus a better approximation, but may require a more points. However, this transformation makes the true integrals simple one dimensional quadratures which can be efficiently computed.</p>
<p>To get the total index from this method, one can calculate the total contribution of the complementary set, i.e. <span class="math">$V_{c_i} = \sum_{j \neq i} V_j$</span> and then</p>
<p class="math">\[
S_{T_i} = 1 - S_{c_i}
\]</p>
<p>Note that this then is a fast measure for the total contribution of variable <span class="math">$i$</span>, including all higher-order nonlinear interactions, all from one-dimensional integrals&#33; &#40;This extension is called extended FAST or eFAST&#41;</p>
<h4>Proof of the Ergodic Theorem</h4>
<p>Look at the map <span class="math">$x_{n+1} = x_n + \alpha (\text{mod} 1)$</span>, where <span class="math">$\alpha$</span> is irrational. This is the irrational rotation map that corresponds to our problem. We wish to prove that in any interval <span class="math">$I$</span>, there is a point of our orbit in this interval.</p>
<p>First let&#39;s prove a useful result: our points get arbitrarily close. Assume that for some finite <span class="math">$\epsilon$</span> that no two points are <span class="math">$\epsilon$</span> apart. This means that we at most have spacings of <span class="math">$\epsilon$</span> between the points, and thus we have at most <span class="math">$\frac{2\pi}{\epsilon}$</span> points &#40;rounded up&#41;. This means our orbit is periodic. This means that there is a <span class="math">$p$</span> such that</p>
<p class="math">\[
x_{n+p} = x_n
\]</p>
<p>which means that <span class="math">$p \alpha = 1$</span> or <span class="math">$p = \frac{1}{\alpha}$</span> which is a contradiction since <span class="math">$\alpha$</span> is irrational.</p>
<p>Thus for every <span class="math">$\epsilon$</span> there are two points which are <span class="math">$\epsilon$</span> apart. Now take any arbitrary <span class="math">$I$</span>. Let <span class="math">$\epsilon < d/2$</span> where <span class="math">$d$</span> is the length of the interval. We have just shown that there are two points <span class="math">$\epsilon$</span> apart, so there is a point that is <span class="math">$x_{n+m}$</span> and <span class="math">$x_{n+k}$</span> which are <span class="math">$<\epsilon$</span> apart. Assuming WLOG <span class="math">$m>k$</span>, this means that <span class="math">$m-k$</span> rotations takes one from <span class="math">$x_{n+k}$</span> to <span class="math">$x_{n+m}$</span>, and so <span class="math">$m-k$</span> rotations is a rotation by <span class="math">$\epsilon$</span>. If we do <span class="math">$\frac{1}{\epsilon}$</span> rounded up rotations, we will then cover the space with intervals of length epsilon, each with one point of the orbit in it. Since <span class="math">$\epsilon < d/2$</span>, one of those intervals is completely encapsolated in <span class="math">$I$</span>, which means there is at least one point in our orbit that is in <span class="math">$I$</span>.</p>
<p>Thus for every interval we have at least one point in our orbit that lies in it, proving that the rotation map with irrational <span class="math">$\alpha$</span> is dense. Note that during the proof we essentially showed as well that if <span class="math">$\alpha$</span> is rational, then the map is periodic based on the denominator of the map in its reduced form.</p>
<h2>A Quick Note on Parallelism</h2>
<p>Very quick note: all of these are hyper parallel since it does the same calculation per parameter or trajectory, and each calculation is long. For quasi-Monte Carlo, after generating &quot;good enough&quot; trajectories, one can evaluate the model at all points in parallel, and then simply do the GSA index measurement. For FAST, one can do each quadrature in parallel.</p>


        <HR/>
        <div class="footer">
          <p>
            Published from <a href="global_sensitivity.jmd">global_sensitivity.jmd</a>
            using <a href="http://github.com/JunoLab/Weave.jl">Weave.jl</a> v0.10.6 on 2020-12-12.
          </p>
        </div>
      </div>
    </div>
  </div>
</BODY>

</HTML>
